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■ Abstract 

I Around 2002, Leonid Gurvits gave a striking randomized algorithm to approximate the 

^ • permanent of an n x n matrix A. The algorithm runs in O (n^/e^) time, and approximates 

U I Per (j4) to within ±e \\A\\" additive error. A major advantage of Gurvits's algorithm is that it 

■ works for arbitrary matrices, not just for nonnegative matrices. This makes it highly relevant 
I to quantum optics, where the permanents of bounded-norm complex matrices play a central 

role. Indeed, the existence of Gurvits's algorithm is why, in their recent work on the hardness 
of quantum optics, Aaronson and Arkhipov (AA) had to talk about sampling problems rather 
than estimation problems. 

In this paper, we improve Gurvits's algorithm in two ways. First, using an idea from 
. quantum optics, we generalize the algorithm so that it yields a better approximation when the 

I matrix A has either repeated rows or repeated columns. Translating back to quantum optics, 

^ ■ this lets us classically estimate the probability of any outcome of an AA-type experiment — even 

an outcome involving multiple photons "bunched" in the same mode — at least as well as that 
probability can be estimated by the experiment itself. (This does not, of course, let us solve the 
I AA sampling problem.) It also yields a general upper bound on the probabilities of "bunched" 

outcomes, which resolves a conjecture of Gurvits and might be of independent physical interest. 
Second, we use e-biased sets to derandomize Gurvits's algorithm, in the special case where 
. the matrix A is nonnegative. More interestingly, we generalize the notion of e-biased sets to 

■ the complex numbers, construct "complex e-biased sets," then use those sets to derandomize 
. even our generalization of Gurvits's algorithm to the multirow/multicolumn case (again for 

nonnegative A). Whether Gurvits's algorithm can be derandomized for general A remains an 
CN ■ outstanding problem. 

■ 1 Introduction 

>< 

■ The permanent of an n x n matrix has played a major role in combinatorics and theoretical computer 
science for decades. Today this function also plays an increasing role in quantum computing, 
because of the fact (pointed out by Caianiello [5] and Troyansky and Tishby [16] among others) 
that the transition amplitudes for n identical, non- interacting bosons are given by n x n permanents. 

Recall that the permanent of an n x n matrix A = (aij), with entries over any field, is defined 
as follows: 



Per(^):= ^ H^^.-W' (1) 
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where 5„ denotes the set of permutations of {1, . . . , n}. A great deal is known about the complexity 
of computing Per (A), given A as input. Most famously, Valiant [T7] showed in 1979 that computing 
Per (A) exactly is #P-complete, even if ^ is a matrix of nonnegative integers (or indeed, if all of 
its entries are either or 1). On the other hand, in a celebrated application of the Markov- 
Chain Monte-Carlo method, Jerrum, Sinclair, and Vigoda [12] gave a randomized algorithm to 
approximate Per (A) to multiplicative error e, in poly (n, 1/e) time, in the special case where A is 
nonnegative. When A is an arbitrary matrix (which could include negative or complex entries), 
it is easy to show that even approximating Per (A) to polynomial multiplicative error remains a 
^P-hard problem; see Aaronson [IJ for example. 

We can, on the other hand, hope to estimate Per {A) to within nontrivial additive error for 
arbitrary A G C"^". And that is exactly what a simple randomized algorithm due to Gurvits 
[8] does: it approximates Per (^) to within additive error ibe||A||" in O (n^/e^) time. Here \\A\\ 
represents the largest singular value of A. The appearance of this linear-algebraic quantity when 
discussing the permanent becomes a little less surprising once one knows the inequality 

|Per(^)|<Pf. (2) 

See Section [2] for Gurvits's algorithm, the proof of its correctness, and the closely-related proof of 
the inequality ([2]). 

In particular, if [/ is a unitary matrix, or a submatrix of a unitary matrix, then \\U\\ < 1, so 
Gurvits's algorithm approximates Per ([/) to within ite additive error in O (n^/e^) time. Equiv- 
alently, the algorithm lets us detect when Per (U) is "anomalously close to 1" (note that, if ?7 is a 
random unitary matrix, then Per (U) will usually be exponentially small, and will accordingly be 
a good additive estimate to Per (U)). This observation makes Gurvits's algorithm highly relevant 
to the field of quantum optics, where subunitary matrices play a central role. 



1.1 Our Results 

In this paper, motivated by the problem of simulating quantum optics on a classical computer, 
we study how far Gurvits's algorithm can be improved. We present two sets of results on this 
question. 

First, we generalize Gurvits's algorithm, as well as the inequality ([2]), to perform better on 
matrices with repeated rows or repeated columns. Let B S c^x'^ be an n x fc complex matrix, and 
let si, . . . , Sfc be positive integers summing to n. Also, let A G be an n x n matrix consisting 

of Si copies of the i*^ column of B, for all i G [k]. Then we give a randomized algorithm that takes 
O {nk/e^^ time, and that estimates Per (A) to within an additive error 

=fc g • —j= \\B\\ . (3) 



Interestingly, we obtain this generalization by interpreting certain formal variables xi, . . . , rr^ arising 
in Gurvits's algorithm as what a physicist would call "bosonic creation operators," and then doing 
what would be natural for creation operators (e.g., replacing each Xi by y/slxi). 
We also prove the general inequality 



si! • • • Sfc! 



Ppr(4)|<- : ^- II (4) 
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Notice that both of these results reduce to Gurvits's in the special case k = n and si = • • • = s„ = 1, 
but in general they can be much better (for example, if we have information about ||i3|| but not 
||74||). We discuss the quantum-optics motivation for these improvements in Section [1.3i 

Second, we show that both Gurvits's algorithm and our generalization of it can be derandomized, 
in the special case where all matrix entries are nonnegative. That is, for nonnegative n x n 
matrices A G M^q", we show that we can approximate Per (A) to within additive error ite • ||^||"', 
deterministically and in poly (n, 1/e) time, li A € ]R"q" is obtained from an nxk matrix B S R>q'^ 
as described above, then we also show that we can approximate Per (A) to within additive error 

±e- \\B\\\ (5) 

deterministically and in poly (n, time. 

To derandomize Gurvits's original algorithm, the idea is simply to use e-biased sets. To deran- 
domize our generalization of Gurvits's algorithm is more interesting: there, we need to generalize 
the notion of e-biased sets to complex roots of unity, then explicitly construct the "complex e-biased 
sets" that we need. 

Let us compare our algorithm to a previous deterministic approximation algorithm for the 
permanent. Gurvits and Samorodnitsky [9] gave a deterministic algorithm for estimating the 
mixed discriminant of a positive semidefinite matrix up to a multiplicative error of e". The 
permanent of a nonnegative matrix is a special case of this. Our algorithm improves over Gurvits 
and Samorodnitsky's [9] when Per (^) is large, i.e., close to ||^||" or to — ^t=^^£t=||i?||"'. 

Of course, for nonnegative A, the algorithm of Jerrum, Sinclair, Vigoda [12] (JSV) gives a better 
approximation to Per (A) than any of the algorithms discussed above (including ours). Crucially, 
though, no one seems to have any idea how to derandomize JSV. One motivation for trying to 
derandomize Gurvits's algorithm is the hope that it might be a stepping-stone toward a derandom- 
ization of the much more complicated JSV algorithm. 



1.2 Quantum Optics 

This paper is about classical algorithms for a classical permanent-estimation problem. However, 
the motivation and even one of the algorithmic ideas came from quantum optics, so a brief discussion 
of the latter might be helpful. For details of quantum optics from a computer-science perspective, 
see for example Aaronson [T] or Aaronson and Arkhipov [2j. 

In quantum (linear) optics, one considers states involving identical, non-interacting photons (or 
other bosonic particles), which can move from one location (or "mode") to another, but are never 
created or destroyed. In more detail, the photons are in a superposition of basis states, each of 
the form |si, . . . , s^), like so: 

1^) = asi,...,sjsi,. . . ,Sk) , (6) 

si,...,Sfc>0 : siH \-Sk=n 

for some complex numbers asi,...,sj. satisfying 

\as,,...,sf = 1. (7) 

si,...,Sfc>0 : siH \-s^:=n 
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Here Si is a nonnegative integer representing the number of photons in the i*'* mode, and si + 
■ ■ ■ + Sfc = n is the total number of photons. To modify the state |$), one apphes a network of 
"beamsphtters" and other optical elements, which induces some k x k unitary transformation U 
acting on the k modes. Via homomorphism, the k x k unitary U gives rise to a ^ C^^n^^) 

unitary ip (U) acting on the n-photon state |<I>). Crucially, each entry of the "large" unitary ip (U) 
is defined in terms of the permanent of a matrix formed from entries of the "small" unitary U. 
The formula is as follows: 

{s^, ...,s,\^{u)\tu..., t,) = («) 

Here Usi^...^SkM,—,tk n x n matrix formed from U by taking Sj copies of the i^^ row of U for 

all i G [fc], and tj copies of the j*^ column of U for all j E [k]. It can be checked [21 [1] — it is not 
obvious! — that is a homomorphism, and that if (U) is unitary for all U. (Indeed, the "reason" 
for the scaling term \/ si\ ■ ■ ■ Sfc'^i' • • • ^fc' is to ensure these properties.) 

Thus, in quantum optics, calculating amplitudes reduces to calculating permanents of the above 
form. By the standard rules of quantum mechanics, the probabilities of measurement outcomes 
can then be obtained by taking the absolute squares of the amplitudes. 

Intuitively, the reason why the permanent arises here is simply that (by assumption) the n 
photons are identical, and therefore we need to sum over all n! possible permutations of photons, 
each of which contributes a term to the final amplitude. Indeed, even in a classical situation, with n 
indistinguishable balls each thrown independently into one of k bins, the probability of a particular 
outcome (for example: 2 balls landing in the first bin, balls landing in the second bin, etc.) could 
be expressed in terms of the permanent of a matrix A of transition probabilities. The difference is 
that, in the classical case, this A would be a nonnegative matrix. And therefore. Per (A) could be 
estimated in randomized polynomial time using the JSV algorithm [12]. In the quantum case, by 
contrast, we want to estimate |Per {A)\'^ for a matrix A of complex numbers — and here is it known 
that multiplicative estimation is already a #P-hard problem [2j. That is why we need to settle 
for additive approximation — or equivalently, for approximation in the special case that |Per(A)|^ 
happens to be "anomalously large" (i.e., non-negligible compared to the general upper bound that 
we prove on |Per (A)|^). 

1.3 Implications of Our Results for Quantum Optics 

With that brief tour of quantum optics out of the way, we can now explain the implications 
for quantum optics of our first set of results (the ones generalizing Gurvits's algorithm to the 
multirow /multicolumn case) . 

Consider the "standard initial state" |1, . . . , 1, 0, . . . , 0), which consists of n identical photons, 
each in a separate mode, and the remaining m — n modes unoccupied. Suppose we pass the n 
photons through a network of beamsplitters, then measure the numbers of photons in each of k 
output modes. Given nonnegative integers si, . . . , Sfc summing to n, let p be the probability that 
exactly si photons are found in the first output mode, exactly S2 are found in the second output 
mode, and so on up to Sk- Then by ([8]), we have 

si! • • • Sfe! 
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where A is the n x n matrix formed by taking Si copies of the i^^ row of the mode-mixing unitary 
U, for all i G [k]. From this, together with the inequality ([1]), we have the upper bound 



P<-ri -ST' (10) 

■^1 "' 



and for the associated amplitude a, 



si! • • • Ski 



V ^1 

As a sample application of p^ . suppose si = n and S2 = • • • = Sfc = 0. Then 

n! 1 

P<— ~— . (12) 

This says that, regardless of what unitary transformation U we apply, we can never get n identical 
photons, initially in n separate modes, to "congregate" into a single one of those modes with more 
than n/e^ probability (the factor of n arising from taking the union bound over the modes). 
Or in other words, for n > 3, there is no direct counterpart to the Hong-Ou-Mandel dip [10], the 
famous effect that causes n = 2 photons initially in separate modes to congregate into the same 
mode with probability 10 

Let us remark that the bound (jlOp is tight. To saturate it, let the n x n unitary U be 
block-diagonal with k blocks, of sizes si x si, . . . , Sfc x Sk- Also, let the i^^ block consist of the 
Sj-dimensional Quantum Fourier Transform, or any other x st unitary matrix whose top row 
equals (l/y^, . . . , 1/y^). Then one can calculate that the probability of observing si photons 
in the first mode of the first block, S2 photons in the first mode of the second block, and so on is 
precisely 

1(1, . . . , 1|99 (?7) 0, . . . , 0, S2, 0, . . . , 0, . . . )|' = '^'"'i . (13) 

^1 " ' ^k 

Here is another implication of our results for quantum optics. Given a description of the unitary 
transformation U applied by the beamsplitter network, the bound ([3]) implies the existence of a 
randomized algorithm, taking O {nk/e^^ time, that estimates the amplitude a = J^^l^,'^'^ \ to within 
an additive error of 




Si! • • • Sfc! 



<^mR5-^ (14) 

V ^1 ■ ■ ■ «fe 

< e, (15) 



I 1 2 

and likewise estimates the probability p = |a| to within an additive error of 

^ I I ; '^J-- "K- \\ -nun ^ Si! ■ ■ ■ Sfe! II oiin f-ic\ 

— i^ll^ll ^^-^^ — ttII^II (16) 
■^1 "' ^k 

< e • 4i 17 

«i • • • ^k 

< e. (18) 




^When n — 2, we get p < 1/2, which corresponds exactly to the Hong-Ou-Mandel dip: the 2 photons have 
probability 1/2 of both being found in the first mode, and probability 1/2 of both being found in the second mode. 
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Observe that the above guarantees match Gurvits's in the special case that ah Sj's are equal to 
1, but they become better than Gurvits's when Y^- \si — 1| is large. In the latter case, (fTO]l says 
that the probability p is exponentially small, but ()16p says that we can nevertheless get a decent 
estimate of p. 



2 Gurvits's Randomized Algorithm for Permanent Estimation 

The starting point for Gurvits's algorithm is the following well-known formula for the permanent 
of an n X n matrix, called Ryser's formula [T5l: 

n 

Per(^)= Yl (-ir+-+""nKi^i + --- + «^-^")- (19) 

x=(xi,...,x„)&{0,l}" i=l 

Ryser's formula leads to an O (2"n^) algorithm for computing the permanent exactly, which can 
be improved to O (2"n) by cleverly iterating through the x's in Gray code order. This remains the 
fastest-known exact algorithm for the permanento 

To present Gurvits's randomized algorithm, we will use a similar formula due to Glynn [7], 
which pulls from the domain { — 1, 1}" rather than {0, 1}". We state Glynn's formula in terms of 
the expectation of a random variable. For a given x G { — 1, 1}", define the Glynn estimator of an 
n X n matrix A as 

n 

Gly^(A) := xi • • • x„ {oi^ixi -\ h ai,nX„) . (20) 

1=1 

Then we have 

Per(^)= E [Gly,(A)]. (21) 

a;6{-l,l}" 

To see why, we just need to expand the product: 

E [Gly^ (A)] = V ai,<^i • • • a„,^„ E [(xi • • • x„) (x^^ ■ ■ ■ x^J] ■ (22) 
xe{-i,i}" ^-^ xe{-i,i}" 

<Ti,...,a„e[n] 

Then, note that E^g|_i [(xi • • • x„)(X(ji • • • Xa„)] is 1 if the map « i— )• cij is a permutation, and 
otherwise. Hence the above sum is simply Per(j4). 

As a special case of a more general algorithm for mixed discriminants, Gurvits [8] gave a 
polynomial-time randomized sampling algorithm for the permanent. His algorithm is simply the 
following: for some m = O (l/e^), first choose n-bit strings xi, . . . ,Xm £ {"l^ 1}" uniformly and 
independently at random. Then compute Gly^_^. (A) for all j € [m], and output 

Gly,^ (^) + ... + Gly,^ (A) ^^^^ 
m 

as the estimate for Per (A) . 



■^Recently Bjorklund _i] discovered a slightly faster algorithm, taking 2" ^(V"/'°s") ^jjjjg^ £q]. computing the 
permanent of a matrix of poly (n)-bit integers. 
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Since each Gly^.. (A) can be computed in O (n^) time, the algorithm clearly takes O (n^/e^) 
time. The algorithm's correctness follows from a single (important) fact. Recall that ||^|| denotes 
the largest singular value of A, or equivalently, 

11.11 ^ „ . \ 

Pll := sup^^—fj^ (24) 

xf^o IfII 

where ||x|| denotes the 2- norm of the vector x. 

Proposition 1. For any n x n complex matrix A, we have |Gly^ (^)| < ||j4||"'. 
Proof. For any x G {1, —1}", we have by the arithmetic-geometric mean inequality, 




|Gly,(A)| < I ^/ ^^=il»^-i^^ + --- + °^."^"lj (25) 

(26) 
(27) 
□ 

To summarize, Gly^. (A) is a random variable that is bounded by WAW^ in absolute value, and 
whose expectation is Per (A) . This implies that 

\FeT {A)\ < \\Af . (28) 

By a standard Chernoff bound, it also implies that we can estimate Per (A) to additive error 
ibe ||j4||", with high probability, by taking the empirical mean of Gly^. (A) for O (l/e^) random x's. 

3 Generalized Gurvits Algorithm 

We now give our generalization of Gurvits's algorithm to the multirow/multicolumn case, as well 
as our generalized upper bound for the permanent. As in Section II. H given nonnegative integers 
si, . . . ,Sh summing to n, we let B S ^"^'^ be an n x /c matrix, and let A G {j^^-x"- i^g iIhq n x n 
matrix in which the i^^ column of B is repeated s-i times. Also, let TZ[j] denote the set of the 
jth j-QQ^g Qf unity, and let X := Tl[si + 1] x • • • x Tl[sk + !]• Then for any x € X, we define the 
generalized Glynn estimator as follows. If a; = (xi, . . . ,Xfc), then let yi = ^/siXi and 

GenGly,(A) := ^ " X v^''' ' ' -V^"' Wvihi + ■■■+ Vkh^k) (29) 

where hij denotes the (i, j) entry of B. We will use this to estimate Per {A), as follows. First we 
will show that 

E [GenGly^ {A)\ = Per {A) . (30) 



Then we will show that 



IGenGly, (^)| < H^f • (31) 
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The sV'" §^ factor in the definition of GenGly^^^ is a normaHzation factor that we need for Lemma 

[2] below to hold. The value of the normalization factor depends on our choice of the scale factors 
when defining the y^'s. Here, we set yi = y/slxi in order to optimize the bound in Lemma [3] below. 



Lemma 2. For the uniform distribution of x over X , the expected value of GenGly^ {A) is 



Proof We have 



E 



E [GenGly, (A)] = Per (A) 



1 I 

-^TTTT^y^'^ • • • Vk'" Ylivi^i,! + ■■■+ ykbi,k 



i=l 



■Si! • • • Ski 



^ " (Ti,...,o-„G[fc] 

Now notice that by symmetry over the roots of unity, 

^Jiy?---yl'')iytr---ytJ]=o, 

unless the product inside happens to evaluate to • • • in which case we have 



E 

xex 



I |2si I 1 2s,, 

|yi| ■■■\yk\ 



" xGX 



\xi\ ■■■\xk\ 



(32) 

(33) 
(34) 

(35) 
(36) 



Therefore, 



E [GenGly, (A)] = s,\ ■ ■ ■ Sk\ 

xGPc 

= Per (A) 



E 



I!''-.'. 



il,...,t„e[fc] : \{i:ti=j}\=Sj Vje[fc]j=i 



(37) 



(38) 



where ([55]) follows from the fact that each product Y[^=i bi,ti appears si! • • • s^! times in equation 

m- ' □ 

Lemma 3. For all x ^ X , 



|GenGly,(A)| < 
Proof. Let y = (yi, . . . , y^). Then 



si! • • • Sk\ 



■^1 ^k 



\B\ 



\y\\2 



|2 I I I |2 

yi\ H v\yk\ 



Hence 



ll^yll < \\B\\\\y\\ = ^\\B\\. 



(39) 



(40) 
(41) 
(42) 

(43) 
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So letting {By) - be the i^^ entry of By, we have 



n (^y). 



i=l 



< 



'n\\B\ 



\B\ 



n 



(44) 



by the arithmetic-geometric mean inequahty. Therefore 



IGenGly, (^)| 



< 



Si Si 
^1 • • • Sfc 

5i! • • • Ski 



Si! • • • Sfc! 



i=l 



\B\ 



(45) 

(46) 
(47) 

□ 

To summarize, GenGly^ (A) is an unbiased estimator for Per (A) that is upper-bounded by 
\\B\\^ in absolute value. This immediately implies that 



I Per (A) I < 



si! • • • Sfc! 



151 



(48) 



It also implies that there exists a randomized algorithm, taking O {nk/e^^ time, that approximates 
Per [A) to within additive error 

+ p- II Rf (49) 



with high probability. Just like in Section [21 that algorithm is simply to choose m = O (l/e^) 
independent random samples xi, . . . , Xm G ^, then output 



GenGly,^ (^) + . . . + GenGly,^ (A) 



m 



(50) 



as the estimate for Per (A) . The correctness of this algorithm follows from a standard Chernoff 
bound. 

As discussed in Section 11.31 in a linear-optics experiment with the standard initial state, the 
above lets us estimate any output amplitude to within additive error i^y^'^i"!^' ^^'^ output 
probability to within additive error ifce "^^i; ■ This is at least as good an approximation as is 
provided by running the experiment itself (and better, if some of the Sj's are large). 



4 Derandomization 

We now discuss the derandomization of Gurvits's algorithm (and its generalization) for nonnegative 
matrices. 
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4.1 Derandomizing Gurvits 



We start by showing how Gurvits's algorithm can be derandomized in the case of computing the 
permanent of a n x n matrix with nonnegative real entries. That is, we give a deterministic 
algorithm that estimates Per (A) to within ite H^H"" additive error with certainty. 
To do this, we need a pseudorandom generator with a specific guarantee. 

Definition 4. We say that a probability distribution D on {0, 1}" is e-biased if for any nonzero 
vector a £ {0, 1}", we have 



E [(-ir-] 



< e. (51) 
Furthermore, we say that a probabilistic algorithm G is e-biased if it outputs an e-biased distribution. 



Note that if G outputs the uniform distribution over {0, 1}", then G is 0-biased. Unfortunately, 
this requires n random bits, and so it takes exponential time to loop over all possible states. We 
will need an e-biased generator with a much smaller seed. Such a generator was constructed by 
Naor and Naor in [14j . who showed the following. 

Theorem 5. There exists an e-biased generator which runs in poly(n, 1/e) time and uses 

0(logn + logl/e) (52) 

bits of randomness. 

This allows us to derandomize Gurvits's algorithm. 

Theorem 6. If D is an e-biased distribution on {0, l}", then for any n x n matrix A with non- 
negative real entries, we have 



Per (A) - E [GlyJ^)] 

x^T> 



<£-\\Ar 



(53) 



where we define the vector x G {—1, 1} by Xi = (—1)*^'. 
Proof. Following equation ([22]) we get 

«n,cr„ E [(Xl • • • Xn){Xai 



E[GlyJ^)]= 

CTi,...,(T„e[n] 



Xr^V 



(54) 



If the map z i—)- cTj is a permutation, then ^er^v {{xi ■ ■ ■ Xn){xa^ ■ ■ ■ Xo-„)] will always be 1. We would 
like to bound the terms for which a does not give a permutation. Given such a o", define a vector 
c G {0, 1}" such that 

Q = l + (|{j:aj=i}|)mod2. (55) 



Clearly, c is nonzero if a does not give a permutation. Then 



E [(xi---x„)(x<^i ■■■x„^)] 



< e 



(56) 
(57) 
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since V is e-biased by assumption. Then, 



Per(A)- E [Gly,(^)] 



o-i,.-.,tT,ig[n] 
i^a^ not a permutation 



< 



< e 



<£ 



E 

t a pc 

E 

jcr: 

E 

■■■,crn€ 

t a per: 



o-i,...,CT„G[n] 
it— fcTj not a permutation 



ai,cri • • • a?! 



CTl,...,o-„G[n] 
i>—^a^ not a permutation 



o-i,...,(T„e[n] 
n 

= ^ • H h a„,j) 

= £-Gly(i_i)(^) 
<e-\\Ar 



(58) 

(59) 
(60) 
(61) 

(62) 

(63) 
(64) 



where for (|59|) we use the fact that the Ojj are nonnegative, and for (|64|) we appeal to Proposition 

m □ 

Corollary 7. There exists an algorithm to deterministically approximate the permanent of annxn 
matrix A with nonnegative entries to within error e ■ \\A\\^ which runs in time poly (n, 1/e). 

Proof. By Theorem [5] there is a polynomial time algorithm which outputs an e-biased distribution 
D on {0,1}". Since the seed for this algorithm is only O (log n + log 1/e) bits, we can compute 
Ee~x) [Glyj.(^)] by looping over all possible seeds. By Theorem [H this suffices as our desired 
approximation. □ 

4.2 Derandomizing the Multicolumn Case 

To derandomize the algorithm of Section [3l we need to generalize the concept of e-bias to arbitrary 
roots of unity. 

Definition 8. Given si, . . . , with n = si+- ■ -+5^, recall that we let X := 7^[si+l] x • • • x7^[sfc+l]. 
Then we say that a probability distribution T) on X is complex-e -biased if for any ei, . . . , e^, with 
ei € {0, . . . , Si\ where not all Cj are 0, we have 



E [x^^ 



< e. 



(65) 



Also, we say that a probabilistic algorithm G is complex-e -biased if it outputs a complex-e -biased 
distribution. 

In Section 14.31 we will show the following: 

Theorem 9. There exists a complex-e -biased generator which runs in poly (n, 1/e) time and uses 
O (log n + log 1 /e) bits of randomness. 
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These complex-e-biased distributions are useful because of the fohowing property, which is 
analogous to Theorem [6l 

Theorem 10. If D is a complex-e -biased distribution on X , then for any n x k matrix B with 
nonnegative real entries, let A be an n x n matrix with Sj copies of the i^^ column of B. Then we 
have 



Per(^)- E [GenGly,(A)] 



< e 



si! • • • Sfc! 



\Br. 



(66) 



Proof. We begin by expanding the product: 

,1 . . . .qJ 



E 




■51 






Si 

Si ■ 




Si!-- 





:yi ■■■Vk 



■Sk 



Jl {yiKi H h ykh,k) 



i=l 
■ br, 



i,<7i • • • on,an E [{yi'^ ■ ■ ■ yk'") (y^i ■■■yaj 

x^T> 



(67) 
(68) 



o-i,...,CT„e[fc] 



-'1,0-1 ' ' ' "n,a„ \ I 



bn.fTr^ A I «^ E [(3^1 - - - Xfi; ^) (^cri * * * -^0 



1 ■■■Sk 



(69) 



Say that a sequence a = (ci, . . . , ct„) is a desired sequence if for all i G {1, . . . ,k}, we have 
|{j : aj = i}\ = Si. Let 

ei = {-Si + \{j : aj = i}|)mod(sfe + 1) . (70) 



Then we have 



(XX - - - X/j *) (Xq-j - - - Xct,j) Xj^ - - - x^ . 



(71) 



Note that = for all i if and only if a is desired. Thus, in the case that a is desired, the expected 
value of ([7T]) is 1. Otherwise, since V is complex-e-biased, we get 



E [(xT^i---x^^'=)(x.,---x,J] 

Xr^T> 



< £. 
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Therefore, 



Per(A)- E [GenGly,(A)] 



Si! • • • Sfc! 



o"! crn&[k] 

a not desired 



'-^^P^ E [(lrT----x^^'=)(x.,---x.J] 



< 



< e 



< e 



si! • • • sk\ 



(7 not desired 



E^[(xT^i---x^^'=)(x,,---x,J] 



Si! • • • Sfc! 



^l"*" ■ ■ ■ Sj!^ (71 CTne[fc] 

a not desired 



Si! • • • Sfc! 



■^1 ■ ■ ■ '^fc cri,...,o-„6[fc] 



CTl -^crn 



Si! • • • Sfc! 



H h ^/Skbi^k) 



e-GenGly(i_...^i)(yl) 



< e • 



si! • • • Ski 



\B\ 



'1 "fc 

where the last inequaUty follows from Lemma EJ 

Theorems [9] and [10] immediately imply the following. 



(72) 
(73) 
(74) 
(75) 
(76) 

(77) 

(78) 
(79) 

□ 



Corollary 11. For any n x k matrix B with nonnegative real entries, let A be an n x n matrix 
with Si copies of the i^'^ column of B. Then we can deterministically approximate Per(yl) to within 
additive error 



si! • • • Sfc! 



• • • sf- 



151 



(80) 



in poly (n, 1/e) time. 



4.3 Constructing Complex-e-Biased Generators 

Our proof of Theorem [9] will be modeled after the proof of Theorem [5] given by Naor and Naor |14j. 
We just need to generalize each step to deal with the domain of complex roots of unity rather than 
the binary domain {1,-1}. 

Recall that our goal is to generate, with as few bits of randomness as possible, complex numbers 
(("i, • • • , Cfc) S A", such that for any exponents ei, . . . , e^, we have 

|E[cr---c]l<^. (81) 
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Figure 1: The ^-strong region of the unit circle in the complex plane. 



In our construction, it will be simpler to adopt an alternate perspective, which we now describe. 

If we let fj be such that Q = e^'''fi/^'i+^\ and let = e2'^*^^/('^^+i), then 

C;^ = (e2-/(^^+i))'^^^=Cf. (82) 

Hence our task is equivalent to outputting a random /c-tuple (/i, . . . , /fc) such that for any choice 
of complex numbers = e^'^*'^i/(*i+^), not all 1, we have 



E 



^{'■■■it] <e. (83) 



The proof of Theorem [5] given by Naor and Naor |14j revolves around the ability to find vectors 
r G {0,1}" such that r ■ x = lmod2 with constant probability for any nonzero x S {0,1}". A 
uniformly random r G {0, 1}" will satisfy r ■ x = 1 mod 2 with probability 1/2, and in fact it is also 
possible to generate such an r using only a small amount of randomness (possibly allowing for a 
weaker constant). 

In our setting, this would translate to generating (/i, . . . , f^) such that Ci^ " " " = ~1 with 
constant probability. Unfortunately, this will not be possible since the do not come from the 
binary domain {1,-1}. Instead, we will have to replace the concept of "being equal to —1" with 
what we call "0-strongness." 

Definition 12. For a unit-norm complex number A, we say that A is 9-strong i/|arg A| > 6, where 
argA is the phase (j) € [— 7r,7r) such that A = e**^. 

Lemma 13. We can generate fi, ■ ■ ■ , fk such that for any ^i, ■ ■ ■ ,6.k, not all 1, we have 

Pr[^f • • • TT/8-strong] > ^ (84) 

using O {log n) bits of randomness. 

The proof of this will require the following fact. 
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Proposition 14. We can generate k random variables {gi, . . . , gk) such that we have gi G {0, . . . , Sj} 

(nearly) uniformly, and c-wise independence between the gi, using 0(c log n) bits of randomness. 



This is well-known. See, e.g., Luby and Wigderson [13]. It is possible to generate k random 
variables ti,...,tk in Fp, where p is a prime of size poly(A:), such that each one is uniformly 
distributed in Fp, and any c of them are independent, using only clog/c random bits. Since our 
desired domains are {0, . . . , Sj}, to get the above result, we have to take the ti mod Si + 1. There 
is a small error since p is not divisible by Sj + 1, but this error is small enough that it does not 
matter for our purposes. 

Another fact we will need was proved by Naor and Naor [1] 



Proposition 15. There is a constant c such that the following holds: Let I and m be integers such 
that m < i < 2m, and let S is a fixed subset of \k] with \S\ = i. If T is a randomly chosen subset 
of [k] such that each element i G [k] is in T with probability 2/m, and the elements are chosen with 
pairwise independence, then 

Pr[|5ur| < c] < ^. (85) 
In particular, they show that we can take c = 7. These facts are enough for the following proof: 

Proof of Lemma [T3[ Suppose that £ > 1 of the are not equal to 1, and let h be the integer such 
that 2^ < £ < 2^~^^. For now, suppose that we know h; we will relax this assumption later. Do 
the following (using Proposition [14]): 

• Generate u = {ui, . . . ,Uk) whose entries are c-wise independent and such that Ui is (nearly) 
uniformly distributed in {0, . . . , Sj}, in the manner stated above. Here, c is the constant from 
Proposition [151 

• Generate w = {wi, . . . , w^) whose entries are pairwise independent and such that vui G {0, 1} 
and FT[wi = 1] = min(l/2'*-^ 1). 



Finally, we define / = (/i, . . . , /fc) by 



° ^!^^ = ° (86) 

Uj it Wi = 1 



Let V := {i : 7^ 1} and W := {i : Wi = 1}. Now we need that with probability at least 1/2, we 
will have 1 < H W\ < c. This follows from Proposition [TJl Now let us examine A := ^(^ ■ ■ ■ ^['^ . 

f fi 

With probability 1/2, we have I < \V r\W\ < c. If this is the case, write A = ^^^'^ " ' ^i'/ 
where 1^ H = {ii, . . . ,id}, and 1 < d < c. Then, the values of fij will be (nearly) uniform in 
{0, . . . , Si^} and independent. Because of the c-independence, for any fixed values of fi^, . . . , fi^, 
we can condition on those choices and still get that /jj is uniform in {0, . . . ,Sij}. Now, we know 

that is an s\^ root of unity other than 1 (but still not necessarily primitive) . Therefore as ^-^^ 
ranges around the unit circle, we can see that 

^h^ (^^^ " " " ^i'd'^ ~ (fixed unit-norm complex number) (87) 

will have a probability of at least 1/4 of being 7r/4-strong (being conservative with constants 
here). Putting it all together, we see that A has a 1/8 probability of being 7r/4-strong, and 
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we have generated A using only O(logn) bits. However, this has ah assumed that we know the 
value of h. Suppose we do not know h. Then, generate O(logn) random bits, and for every 
/i G {0, . . . , [log2 k\} use these bits to construct A*^'^) using the above method. Let /lactuai be the 
value such that 2'^=^'='"^' < (■ < 2'*=^'^*"='' . Now, randomly and indepently choose bits • • • ) ^logjfc 

and let A := (A'-'^^)''''. There is a 1/8 chance that A'-'^'^'^'"^'^ is 7r/4-strong. If this happens to be 
the case, then write 

7^ actual 

If we fix the values of for h ^ /lactuah then the factor on the right becomes fixed, call it A. If 
A is 7r/8-strong, then A will be vr/S-strong if hh^^^^^^ =0. If A is not vr/S-strong, then A will be 
vr/S-strong if ^hactuai ~ ^ (here, using the fact that A^'*''^*"''') is 7r/4-strong) . In either case, there is 
a 1/2 chance that A is vr/S-strong. Combining with the 1/8 probabililty from earlier, we get that 
in total, there is at least a 1/16 chance that A is vr/S-strong. □ 

However, we need a better guarantee than just constant probability given by Lemma [131 In 
particular, we will need to construct a set of A which has high probability (at least 1 — e/2) of 
containing many vr/S-strong elements. We could sample 0(log(l/e)) samples independently, but 
that would require too many random bits. 

To reduce the number of random bits, we use deterministic amplification. The following result 
was shown independently by Cohen and Wigderson [6] and by Impagliazzo and Zuckerman [11], 
using ideas from Ajtai, Komlos, and Szemeredi [3]. The original motivation for these results was to 
amplify error probabilities in BPP algorithms using few bits of randomness. Here is a formulation 
which is useful to us: 

Theorem 16. Let T = {0, 1}^. Then with r + 0{i) bits of randomness we can efficiently generate 
a set yi, . . . ,y£ €z T such that for any S '^T with \S\ > 2 |7~| /3, we have 

Fr[yi £ S for at least i/2 values of i] > 1 - 2'^^^^ (S9) 

This gives the following lemma: 

Lemma 17. There are constants p,q > such that using O (log n + £) bits of randomness, we can 
generate Ai, . . . , A^ such that with probability 1 — 2"^^ , there are at least p£ values of i such that Aj 
is ir/S-strong. 

Proof. We start with a small amount of initial amplification: let t be the constant such that if we 
generate Xi, . . . ,Xt independently as in Lemma [T3| then with probability at least 2/3, one of the A, 
is vr/S-strong, and let r be the number of random bits needed to generate these Ai, . . . , A^. If we 
want to generate Ai, . . . , A^ for some i, then generate groups A^^i, . . . , Xj^t for j G [£/t] according to 
Theorem [T6j Then with probability at least 1 — 2"^^^/*) = 1 — 2~^^^\ at least half of the groups 
will have at least one vr/S-strong member, i.e., at least a l/2t fraction of the A's will be vr/S-strong. 
The total number of random bits used is r + O (^) = O (log n + £). □ 

It is worth noting that Lemma [17] is a slight difference between our construction of complex-e- 
biased distributions and the construction of e-biased distributions by Naor and Naor [Hj. Naor 
and Naor constructed many values, needing only one to be —1, while we construct many values 
needing many of them to be vr/S-strong, for reasons evident in the proof of Theorem [9] Luckily, this 
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was not a problem, since as we saw, the same class of deterministic amplification results applied, 
allowing us to prove Lemma [T71 



Proof of Theorem [3 Generate Ai , . . . , as in Lemma [T71 choosing 

'logi/2(e/2) log^(e/2) 



max 



P 



where /3 := 



1 + 6" 



{0,1}. Then let 



Randomly and independently choose di 



(90) 

,di, each uniformly from 
(91) 



With probability at least 1 — 2 i^ > 1 — e/2, at least p£ of the Xj will be 7r/8-strong. Suppose that 
this is the case. Then 



E [fi] 



=n 

< /3P^ 

< e/2. 



l + Xi 



(92) 

(93) 
(94) 



Thus, this shows that we get that with probability at least 1 — e/2, the norm of the expected value 
of fi is at most e/2. Hence, the total expected value of n, over the random choices of Ai, ■ ■ ■ ,Xe 
and di , . . . , is at most 

(l-|)5 + f(l)<e. (95) 

Now, since we have to choose i values dj, each in {0,1}, this takes O (l) = 0(logl/e) bits of 
randomness. Generating the Ai, . . . , A^ takes O (log n + £) = O (log n + log 1/e) bits of randomness. 
Therefore, the entire procedure then takes O (log n + log 1 /e) bits of randomness. □ 



5 Open Problems 

There are many interesting open problems about Gurvits's algorithm, and the classical simulability 
of quantum optics more generally. Firstly, can we improve the error bound even further in the 
case that the matrix A has both repeated rows and repeated columns? More specifically, can we 
estimate any linear-optical amplitude 

{s^,...,s,\^ {u) \h,...,t,) = ^l^;-;;;;;::';;! (96) 

to ±1/ poly (n) additive error (or better) in polynomial time? This is related to a question raised by 
Aaronson and Arkhipov [2j, of whether quantum optics can be used to solve any decision problems 
that are classically intractable (rather than just sampling and search problems). 

An even more obvious problem is to generalize our derandomization of Gurvits's algorithm so 
that it works for arbitrary matrices, not only matrices with nonnegative entries. Unfortunately, 
our current technique appears unable to handle arbitrary matrices, as it relies too heavily on the 
expansions of Gly(]^_ (A) and GenGly^^^ (A) having no negative terms. 



17 



An alternative approach to derandomizing Gurvits's algorithm would be to prove a general 
"structure theorem," showing that if J7 is a unitary or subunitary matrix, then Per ([/) can only be 
non-negligibly large if U has some special form: for example, if U is close (in a suitable sense) to 
a product of permutation and diagonal matrices. A deterministic algorithm to estimate Per [U) 
could then follow, by simply checking whether the structure is present or not. We conjecture that 
such a structure theorem holds, but do not even know how to formulate the conjecture precisely. 

Even in the realm of nonnegative entries, there is much room for improvement. Among other 
things, it would be extremely interesting to derandomize a Jerrum-Sinclair-Vigoda-type algorithm 
|12j , enabling us deterministically to approximate the permanent of nonnegative matrices to within 
small multiplicative error. 
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